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ABSTRACT 

Aims. This work aims at presenting the first two-dimensional models of an isolated rapidly rotating star that include the derivation of 
the differential rotation and meridional circulation in a self-consistent way. 

Methods. We use spectral methods in multidomains together with a Newton algorithm to determine the steady state solutions including 
differential rotation and meridional circulation for an isolated non-magnetic rapidly rotating early-type star. In particular we devise an 
asymptotic method for small Ekman numbers (small viscosities) that removes the Ekman boundary layer and lifts the degeneracy of 
the inviscid baroclinic solutions. 

Results. For the first time, realistic two-dimensional models of fast rotating stars are computed with the actual baroclinic flows that 
predict the differential rotation and the meridional circulation for intermediate-mass and massive stars. These models nicely compare 
with available data of some nearby fast rotating early-type stars like Ras Alhague (a Oph), Regulus (a Leo) and Vega (a Lyr). It is 
shown that baroclinicity drives a differential rotation with a slow pole and a fast equator and a fast core and a slow envelope. The 
differential rotation is found to increase with mass, with evolution (here measured by the hydrogen mass fraction in the core) and with 
metallicity. The core-envelope interface is found to be a place of strong shear where mixing will be efficient. 

Conclusions. Two-dimensional models offer a new view of fast rotating stars especially of their differential rotation, which turns 
out to be strong at the core-envelope interface. They also offer more accurate models for the interpretation of interferometric and 
spectroscopic data of early-type stars. 
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1. Introduction 

One of the main challenges of nowadays stellar physics is the un- 
derstanding of the global dynamics of stars. Indeed, the new de- 
tailed observations of st ars that report either on magnetic fields 
dPetit etalj|2009ll2010l) or on abund ances patterns in clusters 
dHunter et al.ll2008l: iBrott et al.ll201 lb both crucially depend on 
the fluid flows that pervade the interior of stars. Thermal convec- 
tion has long been related to the generation of magnetic fields 
but the dynamics of convectively stable radiative regions has 
now been considered for at least twenty years as the cornerstone 
of the understanding of surface abundances in stars. Moreover, 
the heart of this large-scale dynamics is the ubiquitous rotation 
of stars. It is now w ell-known that ro tation is a key parameter 
for stellar dynamos (iPetit et al.l I2008I) but it is also crucial to 
the macroscopic transport through stably stratified radiative re- 
gions. There, the baroclinic (or other) flows together with some 
small-scale turbulence determine the so-called rotational mix- 
ing that is now a key pr ocess to interpret stellar abundances 
dMaeder & Mevnetll2000l) . Hence, progress in the understand- 
ing of stars and in the interpretation of the new data sets, which 
come from spectropolarimetry (magnetic fields), spectroscopic 
surveys, interferometry, seismology, require a precise account of 
the distribution and evolution of angular momentum in stars. 

The most common implementation of rotation and rota- 
tional mixin g in stellar evolution codes is the prescription of 
Zahn] ( 1992I) and its subseq uent improvements dTalon & Zahnl 
1997t Maeder & Zahnlll998l) . This modelling designed for one- 
dimensional codes relies on the simplification of differential ro- 
tation which is purely radial (shellular) and of meridional cir- 
culation that is incorporated as an effective diffusion thanks to 



some nat ural hypothe sis on the turbulence in radiative zones of 
stars (see lZahnlll992l) . Although this model has been developed 
for slowly rotating stars, it has been wi dely used, with mor e 
or less success for any rotating star (e.g. lEkstrom et al1l2012l) . 
Clearly, we are missing a more detailed view of the large-scale 
and long-time-scale dynamics of stars, in the first place to put 
bounds on the validity of ID-models. 

A closer view at the global stellar (fluid) dynamics clearly re- 
quires laying aside the spherically symmetric models. Modelling 
dynamics needs more than one dimension, thus a general use of 
two-dimensional models is the obvious next step in stellar mod- 
elling. 

However, jumping into the multidimensional modelling of 
stars is not a sinecure. Attempts have followed one another since 
the first poly tropes of I James! (1 19641) . (see also lRieutordll2006bl 
for a review of the historical landmarks of 2D-models), but 
presently no code can represent stellar evolution in two (or more) 
space dimensions including self-consistently the large-scale dy- 
namics. Most of the present models impose an ad hoc law of 
internal rotation. The reason for that matter of fact is that com- 
puting the baroclinic flows (or others) together with the stellar 
structure is difficult numerically because equations are stiff. 

T hese foreseeable diffi culties motivated the dedicated st ud- 
ies of|Rieutordld2006al) and lEspinosa Lara & Rieutordld2007|) . In 
iRieutordl d2006al) the specific features of baroclinic flows were 
investigated through a simpli fied model using the Boussinesq 
approximation. The next step dEspinosa Lara & Rieutordll2007l) 
has been the computation of a completely radiative star enclosed 
in a spherical container. With this configuration we could com- 
pute self-consistently the flows and the hydrostatic structure of 
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the interior of a fast rotating star. Although this computation al- 
ready gave a good idea of the differential rotation and the merid- 
ional circulation that may be expected in radiative zones of rotat- 
ing stars, the spherical shape of the container enclosing the star 
is not appropriate to give the surface characteristics of such stars, 
which are important from the point of view of observations. In 
addition the computation remained limited to the interior of the 
star since polar pressure was always more than 1CT 5 times the 
central one (realistic values are of order 1CT 12 or less). The nat- 
ural next step that we present here therefore demanded a serious 
effort to overcome the simplification of the previous models. We 
required that the geometry of the coordinates spouses the sur- 
face of the star so that boundary conditions could be correctly 
implemented and that solutions extend up to the true surface of 
the star as given by an optical depth condition. 

While heading to that goal, it turned out that the computa- 
tion of velocity fields was more difficult than ever essentially 
because of the huge variations of density along the radius of the 
star (typically more than eight orders of magnitude). 

The present paper aims at reporting these latest advances in 
the building of the two-dimensional models of rapidly rotating 
stars. In the following section we shall recap our strategy for 
constructing these models and report on the test that were made 
to validate them. In section 3, we discuss the problem of the 
velocity fields and present the boundary layer analysis neces- 
sary for the use of asymptotic solutions. The following section 
presents some preliminary results including a comparison with 
observational data and a first survey of the dependence of dif- 
ferential rotation with global parameters of intermediate-mass 
stars. Some conclusions end the paper. 



There, one recognize Poisson's equation {(p is the gravitational 
potential, p is the density and G is the gravitation constant), the 
equation of entropy S (T is the temperature, v the velocity, F the 
heat flux and e, the nuclear heat sources), the momentum equa- 
tion here written in an inertial frame (P is the pressure and F v 
the viscous force) and lastly the equation of mass conservation. 

These differential equations need to be completed by consti- 
tutive laws. We write the expression of the heat flux as 



F = -Xr^T - ^-TVS 



(2) 



where Xr is the radiative conductivity and Xtmb a turbulent heat 
conductivity. Km =H/M where fi is the ideal gas constant and 
M the mean molecular mass of the fluid. Here, we hypothesize 
that turbulent convection is equivalent to the diffusion of entropy. 
As for the viscous force, the similar level of modelling turbu- 
lence leads us to the use of a turbulent viscosity. We therefore 
assume that 



Av + - V (V ■ v) + 2 (V In fi ■ V) v 

2 

+ V In p x ( V x v) - - ( V ■ v) V In p 



(3) 



where p. is the dynamical viscosity of the gas. 

The equations are also completed by the relations 

[ P = P(p,T) 
\k = k(p, T) 
3 e.(p, T) 



(4) 



2. A recap of the construction of 2D-models 

2.1. Hypothesis 



As in lEspinosa Lara & Rieutord! d2007l) . we restrict our models 
to those describing isolated rotating stars that are in a steady 
state. This idealization, which neglects chemical evolution as 
well as any mass-loss, is of course a first approach of the much 
more complex reality. In addition, we consider the star strictly 
axisymmetric. Breaking of this symmetry comes from turbu- 
lence and magnetic fields. Turbulence is not a problem in prin- 
ciple, since we are interested in long term averages. Effects of 
turbulence should be quasi-steady and axisymmetric (although 
their correct modelling is still a major issue). Magnetic fields 
on the contrary cannot be ignored a priori. Non-axisymmetric 
steady magnetic fields are known to exist in some rotating stars 
(i.e. roAp stars for instance). However, as progress from simple 
to complicate commands it, we first ignore magnetic fields. This 
may not be so bad for some A stars like Vega whose surface 
magnetic fields are rather weak dPetit et al.ll2010l) . 

2.2. Equations 

The partial differential equations that determine the steady state 
of a lonely rotating star are the following: 



A(f> = AnGp 

pTv VS = -V ■ F + £„ 
pv ■ Vv = -VP - pVcp + F r 
V ■ (pv) = 0. 



(1) 



which respectively give the equation of state, the opacities and 
the nuclear heat generation. We recall that in radiative diffusive 
equilibrium, heat conductivity is related to opacity by 



Xr 



\6o-T i 
3i<p 



(5) 



where cr is the Stefan-Boltzmann constant. For the energy gen- 
eration we use an analytical formula that represents hydrogen 
burning either in CNO or pp-chains, namely: 



s t (p, T, X, Z) = e (X, Z)p 2 T- 2/3 exp (A/T l/i ) 



(6) 



as in lEspinosa Lara & Rieutord! (|2007[). It is completed by the 
use of OPAL tables for the computation of the opacity and the 
derivation of the density from the e quation of state (X - 0.7 and 
Z = 0.02 with solar composition of lGrevesse & Noels|[l993l 

2.3. Boundary conditions 

The foregoing partial differential equations require boundary 
conditions for the solutions to be fully determined. These ap- 
ply at the star centre and surface. The conditions at the centre 
are that the solutions should simply be regular. Surface bound- 
ary conditions are more involved. 

The first outer boundary condition is on the gravitational po- 
tential, which should match that of a field in the vacuum, van- 
ishing at infinity. A way of solving this difficulty is the use of 
the analytical solution of Poisson's equation using the Green 
function. T his is the so-called Self-C onsistent Field method de- 
veloped by lOstriker & Mark! d 19681) and used in 2D models of 
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Ijackson et al.l (120051) . We prefer another technique that consists 
in computing this field in a surrounding empty domain bounded 
by a sphere where simple boundary conditions can be applied 
(see iRieutord & Espinosa Laraf2012l) . 

The next boundary conditions are those to be imposed on the 
velocity field. For an isolated star it is natural to use stress-free 
conditions, namely that 

v • n = and ([cr]n) A n — 

where [<x] is the stress tensor. The fluid is assumed to not flow 
outside the star, and to feel no horizontal stress (n is the outer 
normal of the star). 

Finally, the temperature field requires its boundary condi- 
tions. We assume that the star radiates locally as a black body 
and impose 



-XnVT = crT 4 



(7) 



Both of the velocity and temperature boundary conditions 
need to be applied at the surface of the star, which is ill-defined: 
it depends on an arbitrary criterion like the value of an optical 
depth (most commonly). 

We define the bounding surface of the star, i.e. the surface 
where boundary conditions are taken, as the isobar where the 
pressure is 



8po\e 
^pole 



(8) 



where t, is a user-defined optical depth . It is usually set to r, = 
2/3 but t s = 1 is also used dMorell 19971) . On this isobar T = T eS 
at the pole only. Hence, on this surface is not met except at 
the pole. In fact, to determine the temperature field, we just need 
to give the temperature on this surface as a function of latitude. 
For that, we assume that the part of the star that rests above this 
bounding surface can be described by a polytrope of index n. 
Thus the pressure verifies 



Hence we find the temperature profile along the bounding sur- 
face 



gpole K{6) 



!/(«+!) 



(11) 



[geS(9) K pole 

Finally, on the bounding surface we just need to set 

T = T b (G) (12) 

which is quite simple to implement. The polytropic index may 
be chosen from a power law fit of the opacity in the tempera- 
ture/density range spanned by the atmosphere. We indeed know 
that if k cc p e T~ m , then n = (m + 3)/(t + 1). This is not a rig- 
orous description of these layers, which are also baroclinic, but 
it should be a reasonable approximation to determine the funda- 
mental parameters of the star. Nevertheless, if required, a more 
realistic atmospheric model can be used to determine !),(#). In 
the following calculations we shall use r s = 1 and n — 3. Finally, 
we also note that the true equatorial radius is slightly larger than 
the equatorial radius of the bounding surface but the actual dif- 
ference is always less than 10~ 3 and will therefore be neglected. 

2.4. Angular momentum condition 

Even with the foregoing boundary conditions, the problem is 
not fully constrained because the total angular momentum is not 
specified. To complete the formulation of the problem we may 
either enforce the total angular momentum L or specify the equa- 
torial velocity of the star V eq . In the first case we may write 



I r sin OpUip dV 
J 00 



Vip (r = R,Q= n/2) = V eq 
in the second case. 



p=P s (l-Z) 



(9) 2 5 Numerical method 



where £ is a non-dimensional variable measuring the height 
above the bounding surface. As this "atmosphere" is polytropic, 
we have 



T = T b Q-Q 



(10) 



The top of the atmosphere is the true surface of the star as defined 
by a given optical depth t s . There the pressure is 



geff 

P = t s — 

K 



which depends on the co-latitude 6. On the pole obviously p 
P s . The height of the atmosphere is given by 





!/("+!) 


geff Kpole 




P S K _ 




gpole K 



l/(«+l) 



However, at the true surface T = T e ff therefore 

l/Cn+l) 



T e ff - Tb 



geff Kpole 



gpole 



We now briefly re cap on the numerical method that we use. We 
refer the reader to IRieutord & Espinosa Laral d2012l) for a more 
detailed account. 

In order to apply (more) easily boundary conditions we use 
coordinates that spouse the spheroidal shape of the star. The stel- 
lar domain is divided in subdomains. Inside such a subdomain 
the relation between spherical coordinates (r, 9, if) and our new 
spheroidal coordinates ((, 0', if') is simply 



<p = <p 



(13) 



in the i-th domain such that r/, < £ < T]j+i. In this expression, 
Ri(0) is the radius of the i-th bounding surface and A,(^), Bj{£) 
are third order polynomials chosen so that, with the choice 
of t he constants a,- the mapping is continuous and deriv able 
(see IRieutord & Espinosa Larall2012l:lBonazzola et al.ll 19981 for 
more details). Note that these coordinates are the natural coor- 
dinates associated with the shape of the star, but unfortunat ely 
they are not orthogonal (see IRieutord & Espinosa Larall2012l) . 

The spatial discretisation of the equations is done by the use 
of spectral methods so as to minimize the size of matrices that 
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Mass SR/R 
3 1(T 3 
7 6 x lfr ; 



6L/L 
3 x 1(T 3 

3 x lfr 2 



Spc/Pc ST C /T C 

5 x irr 3 8 x icr 3 
5 x irr 2 3 x irr 4 



Table 1. Comparison of the results between the one-dimensional 
of the ESTER code and the one-dimensional code TGEC. 



Req (Rq) 

L(L G ) 

Teq (K) 

Tpo, (K) 
V eq (km/s) 



Deuoree et al. (2012) ESTER 

3.006 2.899 

33.45 33.00 

7735 7798 

9135 9229 

236 238 



Tabl e 2. Compariso n of the results between the models of 
iDeupree et al.l d2012l) and ours for a Ophiuchi. The assumed pa- 
rameters are the mass, M=2.25 M , the chemical composition 
(Z core = 0.25, X env = 0.7, Z = 0.02) and the ratio R poi /R eci = 
0.838. The ESTER model has been computed with 32 spherical 
harmonics and 330 radial grid points spread on 8 spectral sub- 
domains. 



arise. The fields are expanded into a series of spherical harmon- 
ics for their horizontal dependence in the 8' coordinate, and the 
radial functions are sampled on the Gauss-Lobatto collocation 
nodes, which are associated with Chebyshev polynomials. 

The nonlinear equation of the stellar structure are solved it- 
eratively. Two methods have been tested: the first is a relaxation 
method also kn o wn as the fixed-point scheme inspired from 
iBonazzola et al.l (Il998h and lEspinosa Lara & Rieutor 3 (120071) . 
while the second is the well-known Newton's method. The fixed- 
point scheme can be easily coded, while Newton's scheme cod- 
ing is rather involved. However, it turns out that the first scheme 
is really efficient in simple configurations (polytropes or slowly 
rotating stars) while Newton's scheme converges in most cases 
rather rapidly (in typically ten iterations when a non-rotating 
model is given as a starting point). We therefore kept Newton's 
method. 



2.6. Calibrations 

The solutions can be tested in various ways. First, the internal ac- 
curacy is easily accessible with spectral methods since the spec- 
tral expansion readily shows if the spatial resolution is appropri- 
ate or not. A slight change in the initial guess of the solutions 
allows us to test the influence of the round-off errors. 

When the iterations have converged, we test the solution 
in two ways: first with the virial theorem, which relates inte- 
gral quantities derived from the momentum equation and, sec- 
ond, with an energy test that also checks integral quantities but 
from the energy equation. Details of these tests may be found in 
lEspinosa Lara & Rieutord! (120071) or iRieutord & Espinosa Lara! 
(120121) . The virial test is usually passed with and accuracy bet- 
ter than a few 10~ 12 , while the energy test may show errors up to 
10~ 5 . This much looser behaviour of the energy side comes from 
the use of the opacity tables, which generate some noise (either 
of numerical or physical origin). 

Finally, we compared the solutions given by our code at zero 
rotation with solutions given by other codes using completely 
different numerical methods. We show in Tab.[T]one such com- 
parison with the one- dimensional code TGEC (Toulouse-Genev a 
Evolution Code, see iHui-Bon-Hoal I2008L iTheado et all 1201 2l) . 



Simil ar results are obtained with the CESAM code (iMorell 
1 19971) . The small discrepancies between our code and the fore- 
going codes comes from the absence of evolution in our code, 
which allows us using an analytic formula for the heating by nu- 
clear reactions, i.e. ((5). We also find simi lar agreement with the 
values published bv lDeupree et al.l d2012l) although with a slight 
(4%) discrepancy on the equatorial radius (see Tab.|2]l. 



3. The velocity field 

3. 1 . General preliminaries 

The velocity field deserves special attention as it is the solu- 
tion of stiff partial differential equations. In an isolated steady 
rotating star flows come from Reynolds stresses and baroclinic 
torques. The baroclinic torques are created by the misalignment 
of pressure and density gradients and have to be balanced by the 
grad ient of angular v elocity, so that the star rotates differentially 
(e.g. lRieutordft2006d) . 

In an inertial frame the velocity field is governed by 



p(v • V)v = -VP - pV<p + F v 
V ■ (pv) = 0. 



(14) 



If we scale the velocity setting v = 2D.qRw and use R as the 
length scale (fio is an angular velocity to be specified later), we 
can rewrite these two equations in their non-dimensional form: 



p(w ■ V)w = -Vp - pV(p + EpT^w) 
V ■ (pw) = 



(15) 



where the pressure and the potential have been scaled appropri- 
ately. We thus introduced the Ekman number 



PoCIqR 2 



(16) 



where po is a typical value of the dynamical viscosity and po is 
the scale of density. 

The axisymmetry of the solution suggests the decomposition 
of the velocity field into two contributions: the differential rota- 
tion Q(r, 8) and the meridional circulation u(r, 8). We thus write 



w(r, 8) = u(r, 8) + £l(r, 8) x r = u + Clsip 



(17) 



Here, (r, 8, <p) are the usual spherical coordinates and s = r sin 8 
is the distance to the rotation axis. 

Following (fTTT i. the momentum equation can be split into an 
azimuthal part 

s u I 2 8Q. \ 

pu VQ. + 2D.p = Efi A£2+ -— +VO- Vln/i (18) 

s \ s os 



and a meridional part 



p(u ■ V)m - ps£l 2 s = -Vp - pVcp + EpT^u) . 



(19) 



Taking the curl of this expression (divided by p) we obtain the 
vorticity equation associated with the meridional component of 
the velocity field 

dQ. 2 „ Vp x Vp p 



Vx(wxu) - s ip 

oz 



+ E-Aco + EMJu) (20) 
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Fig. 1. Top: Profiles of the radiative (black) and turbulent (red) 
viscosity for 3 stellar models of ZAMS stars of 3,7 and 10 M . 
Bottom: the associated profiles of the Ekman number v/2QR 2 
based on the radius of the star and a rotation period of 12 hrs. 



where Ri c is the critical Richardson number. Assuming that 
Ri e ~ 1 1 A, we show in Fig. Q] the values of the turbulent Ekman 
number, with the same scalings as above, computed from the dif- 
ferential rotation of a 2D-model. Although a few orders of mag- 
nitude larger, the turbulent Ekman number of fast rotating stars 
is still very small compared to unity (always less than 10" 7 ). 

3.3. Inviscid solution and the occurence of a boundary layer 

The extremely small values of the stellar Ekman numbers lead 
us to look for solutions of the velocity equations valid in the 
asymptotic regime E — > 0. In this regime, the viscous force is 
not able to maintain a solid body rotation. It therefore carries 
angular momentum. In a steady state, this flux is balanced by the 
meridional flow according to ( l22l . This latter expression gives 
the order of magnitude of this meridional flow. Thus, we have 



Ml ~ O(E) 



(23) 



Using this result, the meridional flow can be decoupled from the 
rest of the equations. Indeed, all the terms depending on u in 
the other components of the momentum equation are 0(E 2 ) and 
can be neglected. In particular, the meridional projection will 
become 



psQrs = Vp + pV0 

and the vorticity equation 



(24) 



with 



2(Vlnp ■ V)h 



+Vln^/ x (V x u) - - (V ■ u) Yin// 



(21) 



where oj = V x u is the vorticity of the meridional circulation. 
Note that only the ^-component of this equation is non-zero. 

Eq. <TT~8T > can be easily transformed into an equation for the 
transport of angular momentum 



V ■ (ps 2 Qu) = EV ■ OurVn) 



(22) 



establishing that, for the fluid to be in a stationary state, the 
transport of angular momentum by the viscous force must be 
balanced by the transport induced by the meridional circulation. 

3.2. Viscosity in stars 

For the stars we are considering that have no convective enve- 
lope, the physical conditions are such that radiative viscosity 
largely dominates over that of collisional origin. In Fig. \T\ we 
plot the values of this viscosity along with the associated Ekman 
number based on the radius of the star and an angular velocity 
corresponding to a rotation period of 12 hours. This figure shows 
that the radiative viscosity leads to Ekman numbers always less 
than 1(T 10 . 

Because of these low values, some turbulence may develop 
as a result of the shear of the differential rotation. IZahnl d 19921) 
proposes that the induced turbulent vertical viscosity reads 



Ri c K 



s dQ. 
N~d7 



dO 2 „ Vp x Vp 

S— = <f ■ 7, 



dz 



(25) 



In principle, ( l25l ) can be used to calculate the rotation profile for 
a given configuration. Unfortunately, without adequate boundary 
conditions, the rotation profile is undetermined by a function of 
s, namely if Q(r) is a solution, Q' 2 (r) = Q 2 (r) + F(s) is also a 
solution. 

This degeneracy is removed by viscosity. Indeed, the free 
surface of the star is assumed to support no stress. Hence the 
angular velocity should verify: 



vn = o 



(26) 



where h is the unit vector normal to the surface. However, this 
condition is incompatible with (l25l l. This is a standard problem 
that is solved using a boundary layer analysis. The main idea is 
to take advantage of the existence of a thin layer in the vicinity 
of the boundary where viscosity is not negligible and where the 
flow adapts its properties to fulfil the boundary conditions. Thus 
doing, we can separate the problem in two domains, an interior 
one where viscosity can be neglected and a thin boundary layer, 
dominated by viscosity. In our case, the thickness of the bound- 
ary layer depends on the Ekman number. It is infinitesimally thin 
when E — > 0, which is used to simplify the equations, as the ver- 
tical gradient of a boundary layer quantity is much larger than 
the horizontal one. This procedure has been successfully used 
to derive barocl inic flows in a Boussinesq model of a star (see 
lRieutordll2006al) . We now adapt it to realistic stellar models. 

3.4. Stretched coordinates 

It is useful to define a new set of stretched coordinates (£, r, <p) to 
represent the fluid inside the boundary layer. The vertical coordi- 
nate £ measures the distance to the surface scaled by a parameter 
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Fig. 2. A schematic view of the local basis adapted to the analy- 
sis of the Ekman layer. 



e that represents the thickness of the boundary layer. The angu- 
lar coordinate t is equal to the colatitude of the nearest point on 
the surface, so that r = 9 when £, = 0. Choosing that £ increases 
outwards, £ is negative inside the star. The azimuthal coordinate 
(f coincides with the corresponding spherical coordinate. 

The unit vectors £ and f are respectively perpendicular and 
parallel to the surface. They can be expressed as 



£ = s sin(r - S) + z cos(t - 6) 
f = s cos(t -5)-z sin(r - 8) 



(27) 



where s and z are the unit vectors associated with the cylindrical 
coordinates; 6 is defined as 



tan 6(t) = 



1 dR 



(28) 



and accounts for the deformation of the surface due to the cen- 
trifugal force. Fig. [2] sketches out this new local frame. The rela- 
tion between , t) and the cylindrical coordinates (s, z) is given 
by 



s = R(t) sin t + e^s ■ i; = R(t) sin t + e£ sin(r - 5) 
z = R(t) cos t + e£z ■ £, = R(t) cos t + e£, cos(r - S) 



(29) 



The coordinates (£, r, <p) form a set of orthogonal coordinates 
with scale factors 



h = 

h T = 



h v = 



dr 

d£ 
dr 

dr_ 



cos 5(t) \ dr 



= R(t) sin t + s£ sin(r - 5) 



(30) 



We are interested in a very thin layer below the surface, thus 
s «: 1 and to leading order 



R(t) 



h T = 

cos 5{f) 

h v = R(t) sin t 



(31) 



3.5. Scale analysis 

We start the analysis of the boundary layer by deriving the or- 
der of magnitude of the dependent variables. As we have noted 
before (Eq. [23), the meridional velocity in the interior inviscid 
domain is u ~ 0(E), while the rest of the variables are supposed 
to be ~ 0(1). Since the flow is steady, mass conservation im- 
poses 



u=0 at £, = 



(32) 



Using the boundary layer coordinates, it reads — £, ■ u — 0. So 
the normal velocity should go from a value 0(E) at the base of 



On. 



the boundary layer to zero. Thus ~ 0(E) since £ is a stretched 
variable. Using the new coordinates (£, t, tp), the continuity equa- 
tion now reads (see appendix fAl 



ia_ 



(pwf ) + 



cos 6 d 
~R~~{h 



(pu T ) 



cos(r - 6) 
R sinr 



-pu T 







(33) 







, dQ, 
and — 

a? 



o 



(34) 



from where we see that u T ~ 0(E/s). 

At the surface, the tangential stress should vanish, according 
to the stress-free boundary condition. Using the boundary layer 
coordinates this condition reduces to 
du T 

Let's go back to the vorticity equation (1201) . If the vertical gradi- 
ent of angular velocity vanishes at the surface, the inviscid vor- 
ticity equation d25l l cannot, in general, be verified. Thus, within 
the boundary layer, the balance of forces requires the viscous 
force whose dominant term E^Ao> should be of order unity. The 
vorticity of the meridional velocity field is 

1 du T 

and 

E—Aco = E-- fiur 
P 



co = V x u 



cos 5 dug \ 



0(E 2 /s 4 ) , 



1 du T „ 



pe 3 <9£ 3 

then the scale of the boundary layer must be 
£= Ve. 



(35) 



(36) 



(37) 



as expected. Normal and meridional velocities are respectively 
0(s 2 ) and (9(e). It is useful to define the scaled boundary layer 
corrections 



U = 



and V = 



U T — T ■ Hj, 



(38) 



where u ml ~ 0(e 2 ) is the interior meridional circulation. 

Back to the meridional component of the momentum equa- 
tion ( fT9l , the viscous force is 0(s 2 ) in the direction of £■ and (9(e) 
in the direction of f and, in the limit e — > 0, the inviscid equation 
is still valid within the boundary layer. Then we suppose that the 
rest of the dependent variables p, p, <p are not perturbed in the 
boundary layer and can be treated as functions of r only, as their 
vertical derivatives ^ are (9(e). 

Q is slightly perturbed within the boundary layer: its vertical 
derivative should change from some finite value to zero, thus we 
set 



<9Q _ /dQ\ 
¥ " \d£k 



+ ep\ 



(39) 



where the vertical variation of /? is not negligible inside the 
boundary layer. 
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s (R s ) s (R e ) 



Fig. 4. Streamlines of the meridional circulation for the inviscid solution with E — > (left) and for the full problem (E = 10" 7 ) 
including the viscous force (right). Solid lines and dashed lines represent counter-clockwise and clockwise circulation respectively. 
Note that in the viscous solution the streamlines are closed inside the boundary layer. 




Fig. 3. Surface rotation (in cycles per day) as a function of colat- 
itude for different values of the Ekman number. The profile ap- 
proaches the inviscid solution (dashed line) as the Ekman num- 
ber decreases. 



3.6. Boundary layer equations 

In short, we have to solve for three variables in the boundary 
layer: the boundary layer corrections of the normal and tangen- 
tial components of the meridional velocity field U and V, and of 
the vertical gradient of rotation that we have called /3. The equa- 
tions describing the flow in the boundary layer are obtained from 
the vorticity equation d20b . the conservation of angular momen- 
tum dTST ) and the continuity equation V • (pu) = 0. At leading 
order in s and after subtracting the interior inviscid equations, 
these equations read (see appendixlAl 



d 3 V 

2Qiflcos(T-5) + v— -5- = 0, 
2£l + s— \Vcos(t-6)-sv£ = 0, 

OS j 01; 

s 8U dV I sdp\ 

z- + s—+\l + -^-\V = 0, 

cos(r - 6) os \ p os J 



(40) 
(41) 
(42) 



where we have defined s = R(j) sin t and v = nip. s is the radial 
cylindrical coordinate at the surface and v is the non-dimensional 
kinematic viscosity. These equations are completed with bound- 
ary conditions. At the surface = 0), the conservation of mass 
and the stress-free boundary conditions become 

= -!?#• "int, ^=0 and /3=Mt) at <f = (43) 

E Ot; 



where fi = -£ ■ (Vfl);, 
the interior, then 



The corrections should also vanish in 



U -> , V ^ and /3 -> when f 
Combining ( |40T > and ( BTb . we obtain 



(44) 



d 4 V 1 



— + -r 4Q Z + S— COS^T - 6)V = , 



ds 



(45) 
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whose solution is 

V = VbGDe* [cos(B^) + V!(S)8in(2*0] , 
with 



(46) 



(47) 



Here, we have used the fact that V should vanish when £ — » -oo. 
Using the boundary condition at the surface d43l i. we see that 
V\ = - 1 and therefore 



y = y (S)e Bf [cos(fi£) - sin(fi£)] . 



(48) 



In the previous expressions we have considered that B is real. 
For that, Q. must satisfy the condition 



8Q? 4Q 2 
-^z~ + — >0, 

os s 

which is equivalent to 



(49) 



(50) 



In other words, the specific angular momentum must increase 
with the distance to the axis of rotation. This is in fact Rayleigh's 
criterion for the centrifugal stability of the differential rotation 
profile. 

Using we get B 



B = B (s)e B f cos(B& , 

and the condition 

Q5cos(t-£) 

Vo(i) = T55 ^ o(i) 

2B i v 



(51) 



(52) 



Finally, using the continuity equation (l42l i we obtain the expres- 
sion for U 

U = -C/o(i)e Bf cos(B^)-C/i(S)^e Bf [cos(B^)-sin(B^)] , (53) 
where 



cos(r-c5) d jspV \ 
U ° = sp dil~B~i 



and 



t/i 



Vq cos(t - 5) dfi 



5 



d.v 



(54) 



(55) 



The foregoing expression of the boundary layer correction 
are those of the classical Ekman layer adapted for a spheroidal 
stress-free surface. As the Ekman layer, these corrections suf- 
fers of the equatorial singularity. This singularity is associated 
with the vanishing of B(s) at equator, which is the inverse of the 
layer's thickness. This is a quantity of order unity except at equa- 
tor where it vanishes (r — > n /2 and 6 — > 0). As was shown by 
iRoberts & Stewartsonl (1 19631) . the boundary layer changes scale 
in this region which extends as 0(E 1 ^ 5 ) in latitude. There, the 
boundary layer thickness is 0(E 2 ^ 5 ), which is slightly larger than 
E 1 ^ 2 . The influence of this singularity tends to zero when the 
Ekman number vanishes. 



3.7. Surface differential rotation 

The normal velocity in the boundary layer ( T53l > should verify the 
boundary condition at the surface, as expressed in d43l , namely 



U(£ = 0) 



-Uo = ~—g ■ "int. ■ 

E 



(56) 



We shall see now that this leads to the appropriate boundary con- 
dition on the Q. profile at the surface. As the fluid flow is axisym- 
metric the meridional velocity may be described by a stream 
function if/, defined as 



pu int = V x (iff$) , 
then 

cos(r - 6) 1 d 

Uo = z -tt (stfr) . 

sp E os 

Using d54b we see that, at the surface 
Mr = R(6),6) = Ep^ . 

D 

Using ( |47| ) and d52i l we calculate ^ 



(57) 



(58) 



(59) 



Vb 



B (20 + ~s§) cos(r - S) * ' ^(S 2 ^) 

Eqs. ( |59l > and d60b lead to the condition for the interior rotation 
profile, now expressed using interior variables only 



Efis 2 i ■ VQ. + ifff ■ V(s 2 Q.) = at the surface . 



(61) 



This boundary condition is to be imposed to the differential ro- 
tation and lifts the indetermination of the inviscid solution of 
the baroclinic flow (1251 1. Indeed these latter solutions are deter- 
mined up to an arbitrary function of the radial cylindrical co- 
ordinate s. This new condition leads to the differential equation 
that determines this arbitrary function. It couples the horizon- 
tal derivative of the angular velocity f • V(s 2 Q) with its vertical 
derivative £■ ■ V£X Thus, solving for this nonlinear boundary con- 
dition together with the equations of dynamics directly leads to 
a self-consistent differential rotation without having to solve the 
Ekman boundary layer. This latter part of the solution may be 
recovered using d48l i, (IBTt and (l53l l. 

Remarkably enough, this boundary condition is smooth at 
equator: the above mentioned equatorial singularity is removed 
with the rest of the Ekman layer. 

3.8. Numerical validation 

To validate the new boundary condition d6Tb to be imposed to the 
inviscid baroclinic flows, we solved numerically the full viscous 
equation in the simplified set-up of a star en closed in a spherical 
box as in lEspinosa Lara & Rieutordl d2007l) . As in this previous 
work we considered a one solar mass of an ideal gas, whose 
opacities are given by Kramer's power laws and which rotates 
at half the break-up angular velocity. The dynamical viscosity of 
the gas is set to a constant p c . The interior is fully radiative and 
the polar pressure is set to 2 x 10 the central pressure. This 
latter condition reduces the density variations and makes the full 
viscous solution more easily computable. 
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Bi Bi B£ 

Fig. 5. Boundary layer corrections for the normal velocity (left), latitudinal velocity (centre) and normal gradient of angular velocity 
(right) for different values of the Ekman number E. The corrections are computed at a colatitude 9 = n/4- and are in normalized 
units. Solid lines show the corrections calculated from the difference between the full solution and the inviscid model. Dashed lines 
show the corrections as predicted by the asymptotic solutions (l48l l. (|5T1 i and (153) . 



Because of the imposed constant dynamical viscosity, den- 
sity variations induce variations of the kinematic viscosity and 
thus of the Ekman number ( fTol i. In Fig. [3] we show the conver- 
gence of the surface differential rotation towards that prescribed 
by (|6TT >. when the (central) Ekman number E c = /j./p c 2QR 2 de- 
creases from 1(T 4 to 1(T 6 . In Fig. g] we show the meridional 
stream function in the asymptotic case and for a viscous solu- 
tion with E = 10~ 7 . We note that the two solutions nicely agree. 

Finally, we computed the boundary layer corrections for the 
three components of the velocity both from the full numerical 
solution and from the asymptotic formulae (148b , ( BTT l and i53i . 
As shown by Fig. [5] we see that, provided that the Ekman num- 
ber is sufficiently low, the two solutions also nicely agree. 



chemical evolution of the core), forces a strong angular veloc- 
ity gradient at this same place. In our models, which neglect 
viscous effect in the interior and any diffusion of elements, this 
velocity gradient is represented by a jump of the angular veloc- 
ity. This jump is controlled by the latitudinal variations of the 
pressure along the CEI. Indeed, the CEI is a surface where the 
Brunt-Vaisala frequency squared changes sign. If there is a den- 
sity jump there, N 2 behaves there as a Dirac and the baroclinic 
torque is not defined. However, since pressure must be contin- 
uous, its latitudinal variations on the CEI just inside and just 
outside the core must be the same. From the projection of (l24l) 
on a tangent vector ej of the CEI we can therefore derive the 
following jump condition 



4. Results 

4.1. The Core-Envelope Interface 

The Core-Envelope Interface (thereafter CEI) is the locus of a 
complex dynamics that is important to understand for its conse- 
quences on the lifetime of the stars and the chemical enrichments 
of their surface layers. 

As a consequence of the baroclinic torque balance d25l ). a 
strong density gradient at the CEI (driven for instance by the 



Pcoie(f) (sQ.l OK (ff)e s -e T -e T - V0 core ) = 

Penv(6») (snl m (9)e s -e T -e T - V0 env ) (62) 

This interface conditions shows that if the density is discontinu- 
ous, so is the angular velocity since V<p is continuous. 

Numerics show that if the core is denser than the envelope 
it rotates faster as illustrated in Fig. [6] We note that this faster 
rotation of the core and the ensuing velocity gradient builds up 
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R (R B ) 

Fig. 6. Angular velocity (in rd/s) as a function of the radial dis- 
tance (in R ) along the rotation axis for the Vega model of Tab. [4] 
The angular velocity discontinuity at the core-envelope interface 
is clearly visible. 



on the thermal time scale of the star. This is indeed the time 
scale that controls the baroclinic modes and therefore the action 
of the baroclinic torque. This time scale is usually smaller than 
that of nuclear evolution and therefore this velocity jump has 
time to build up when the core gets denser as nuclear evolution 
proceeds. 

The foregoing shear is obviously smoothed by diffusion pro- 
cesses in particular by viscosity. However, this smoothing is 
not simple. The balance of forces requires the occurence of a 
Stewartson layer that develops on the tangent cylinder of the 
core-envelope interface as illustrated in Fig. [7] Such a layer is 
triggered by any discontinuity tha t appears on this in terface. As 
shown in the original analysis of IStewartsonl (1 19661) . this inter- 
nal layer is composed of nested layers whose thicknesses scale 
with some fractional power of the viscosity (E 1 ^ 4 , E 2 ^ and E 1 ^ 
are the main non-dimensional scales). Such scalings are not in- 
cluded in our solution, which only retains the balance between 
viscous force and meridional advection at the global scale. 

In Fig. 13 we clearly see the meridional flows induced by the 
core-envelope interface. Their spectral expansion on spherical 
harmonics is not well-converged as a consequence of the hidden 
discontinuities. The foregoing remarks aim at underlining that 
such flows should be considered as not fully computed. This is 
why we shall not comment on the transport they may achieve 
although this is one of the key part of rotational mixing. The 
core-envelope shear layer is clearly reminiscent of the tachocline 
of low mass stars, although the mechanisms at work are different. 
As for this layer, we should expect some dynamo effect driven 
by the shear. Further work is clearly needed. 

Fortunately, the differential rotation is much better com- 
puted and does not suffer from the approximate treatment of the 
Stewartson layer. This is most probably because the influence of 
this layer vanishes as the viscosity vanishes. 

4.2. Comparison with observational data 

The most important test of the foregoing models naturally comes 
from a comparison with observational data. The recent pro- 
gresses of optical and infrared interferometry lead to the deriva- 




Fig. 7. Meridional streamlines from the model of a Leo. The 
vertical solid line materializes the tangent cylinder associated 
with the core-envelope interface. 

Tabl e 3. Compari son between observationally derived 
dMonnier et al.l 1201 Ol) parameters of the star a Oph and our 
model. The mass fraction of hydrogen in the core is related to 
the age of the star. P eq and P po i are the rotation period in days at 
equator and pole as predicted by the model. A solar composition 
(X=0.7, Z=0.02) is assumed for the envelope. 





Observations 


Model 


Mass (Mq) 


o 4+0.23 
*• -0.37 


2.22 


Req (Rq) 


2.858±0.015 


2.865 


Rpol (Ro) 


2.388±0.013 


2.385 


Tec, (K) 


7570±124 


7674 


T P „,(K) 


9384±154 


9236 


L(L Q ) 


31.3+1 


31.1 


V eq (km/s) 


240±12 


242 


P eq (days) 




0.598 


Ppoi (days) 




0.616 


Xcore/Xenv. 




0.37 



tion of new precise measurements of the fundamental parameters 
of some early-type stars of intermediate mass. 

We have selected three such stars, namely a Lyr (M ~ 2.2 
M Q ), a Oph (M ~ 2.2 M ) and a Leo (M ~ 4. 1 M H ). a Oph 
(RasAlhague) was also successfully modeled bv lDeupree et al.l 
(120121) . 

Data from these three st ars comes from the interferometric 
observations with CHARA dAufdenberg et al1l2006t IZhao et aTl 
120091: IChe et alJl201 it IMonnier et alJl2012h . In addition, a Oph 
is the primary of binary system whose orbit constrains its m ass 
to the interval 2.03 < M/M < 2.63 dHinklev et alJl201ll) . It 
has also been recently identified as a 5-Scuti star displaying 
57 frequencies ( IMonnier et alj|20loh . All these data make Ras 
Alhague an interesting test bed for models of rapidly rotating 
stars. Vega being a photo metric standard, it has been studied in 
detail for many years (e.g. Aufdenberg et al. 2006; Takeda et al. 



120081 : IMonnier et all 12012 . and references therein). Regulus is 
also an interesting case because it may owe its fast rotation to 
the swallowing of the enve lope of its now white dwarf compan- 
ion dRappaport et alJl2009l) . 
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Fig. 8. Convectively unstable layers in the outer layers of Vega's model. Left: a global view showing the regions where N 2 < 
These are limited by the solid lines. The dashed line marks the minimum of A^ 2 . The layers associated with helium ionization region 
clearly appear. Left: Closer view of the surface layer. 



Table 4. Same as in Tab.[3]but for Vega (a Lyr). O bservationnaly 
derive d values are from the concordance model of lMonnier et alJ 
J20l2l) . A sub-solar met allicity of X=0.75 46, Y=0.2361 and 
Z=0.0093 is adopted after lYoon etall d2008l) . 





Observations 


Model 


Mass (Mq) 


9 1 c+0.10 

A-1 -0.15 


2.374 


Rcq (Ro) 


2.726±0.006 


2.726 


Rpoi (Ro) 


2.418±0.012 


2.418 


T eq (K) 


8910+130 


8973 


Tpoi(K) 


10070+90 


10070 


L(L Q ) 


47.2±2 


48.0 


V eq (km/s) 


197±23 


205 


P eq (days) 




0.672 


Ppoi (days) 




0.697 






0.271 



Table 5. Same as in Tab.|3]but for a L eo. Observationnally de- 
rived values are from (Che et al. 201 1). 





Observations 


Model 


M(M ) 


4.15±0.06 


4.10 


Rcq (Ra) 


4.21±0.07 


4.24 


Rpoi (Ro) 


3.22±0.05 


3.23 


T cq (K) 


11010±520 


11175 


Tp„, (K) 


14520±690 


14567 


L(L Q ) 


341+27 


351 


Veq (km/s) 


336±24 


335 


P C q (days) 




0.641 


Ppoi (days) 




0.658 


Xcore/Xenv. 




0.5 



In Tab.[3j0and[5]we confront our models to these data. Three 
parameters have been adjusted: the mass, the angular velocity 
and the mass fraction of hydrogen in the core. This latter pa- 
rameter allows us to take into account the hydrogen depletion in 
the core due to time evolution. This is our proxy for time evo- 
lution on the main sequence. A solar metallicity is used for the 
envelope (X=0.7 and Z=0.02), ex cept for Vega whe re we use 
Z=0.0093 and X=0.7546 following lYobn et all(l2008l) . 



The comparison shown by these tables is quite encouraging. 
We shall not discuss the matching between data and models in 
more detail. Indeed, the observed quantities are model depen- 
dent: actually, they have been derived using a Roche model with 
solid body rotation and ad hoc gravity darkening laws. A full 
discussion, which should include the spectral energy distribu- 
tion and the oscillation spectrum when relevant, would require a 
dedicated work and is therefore beyond the scope of the present 
paper. 

4.3. Surface convection 

While the most important convective region for stars with a 
mass above two solar masses is the core, some residual con- 
vection may affect the surface layers. In Fig. [8j we display the 
regions where the squared Brunt-Vaisala frequency is negative 
for a Vega-like model. It shows that a thin convective layer af- 
fects the surface due to the hydrogen ionization peak. Two other 
layers associated with the first and second ionization of helium 
are also shown. As an effect of rotation, these layers thicken and 
deepen in the equatorial region. In the more massive star a Leo, 
the unstably stratified layers are below the surface in the polar 
regions, but the equatorial temperature drop lets the outermost 
layer almost touch the surface (see Pig . |9j> ■ These regions, if they 
actually develop some thermal convection may participate to the 
line broadening with micro- or macroturbulence. 

4.4. Differential rotation of intermediate-mass stars 

The essential progress that has been accomplished in these 2D- 
models with respect to previous attempts is that differential ro- 
tation is no longer arbitrary. It is the solution of fluid dynamics 
equations and is the response to the actual baroclinic torque that 
exists in all rotating stars. When stars rotate rapidly and do not 
lose angular momentum, this torque is the main driver of the 
differential rotation. Reynolds stresses from the convective core 
are also possibly influential, but presently their functional form 
is unknown. 

In order to have an idea of the resulting differential rota- 
tion in intermediate-mass stars, we explore in a rather system- 
atic manner the dependence of differential rotation with global 
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Fig. 9. Layers with unstable stratification in the upper envelope 
of the model for a Leo. 




' ", ''»' '<?' 




Fig. 10. Differential rotation of Regulus. Dashed isocontours and 
dark background represent slower regions, while clear or white 
background and solid isocontours show faster regions. 



quantities that may influence it. We therefore examine the de- 
pendence with mass, metallicity and angular velocity. 

We first note that in all cases investigated, this differential 
rotation is stable with respect to the centrifugal instability: the z- 
component of the specific angular momentum always increases 
with the distance to the rotation axis, except for small jump at 
the CEI when the core is evolved. This is a consequence of the 
weak differential rotation generated by the baroclinic torques as 
shown below. 

A first general view of the differential rotation is given in 
Fig. [10] for the model of Regulus. As emphasized by the surface 
latitudinal and radial variations on a ZAMS (homogeneous) 3 so- 
lar mass star model, the differential rotation never exceeds 7 % 
horizontally and 25 % radially (see Fig, [m and [TBI. As shown in 
Fig-HJJthe differential rotation is not necessarily a monotonous 
function of the latitude and at fast rotation latitudes at ±40° 
rotate slightly more rapidly ( 2 %) than the equator. This is a t 
odds with our previous result (lEspinosa Lara & Rieutordll2007h . 
The difference likely comes from the unrealistic boundary con- 
ditions used in lEspinosa Lara & Rieutord! ((2007) where the star 



was confined in a spherical container. The non-monotonous be- 
haviour is present only in the outer part of the envelope. It is 
related to the change of the temperature variation with latitude 
on an isobar, since these variation s control the baroclinic torque 
(lEspinosa Lara & Rieutordll2007l) . Indeed as shown in Fig.[T2"lin 
the interior the temperature decreases from pole to equator while 
near the surface it increases. Such a change affects the baroclinic 
torque and the ensuing velocity field. The reason of this change 
is not clear. 

The radial profiles displayed in Fig.[l3]show the almost rigid 
rotation of the core (as there is no baroclinic torque there). As 
noted above the radial variations are stronger than the latitudinal 
ones, but they are essentially concentrated near the core (namely 
between r c and 2r c ). Since the latitudinal differential rotation 
is much less, we note that in this neighbourhood of the core, 
the differential rotation is essentially shellular (see also Fig.lTOl. 
This remark prompted us to compare the radial differential with 
that obtained from the previous ID -models including rotation 
along the approach of Zahnl (1 19921). Such works like those of 
iDenissenkov et al.l d 19991) and lMevnet & Maeded d2000t) showed 
that an initial solid body rotation rapidly relaxes towards a quasi- 
steady shellular differential rotation like the one displayed in 
Fig. [13] Our model shows that this shellular flow is just the dif- 
ferential rotation driven by the baroclinic torque of the config- 
uration and does not depend on a prescription for the viscosity. 
As a consequence the ID and 2D model are in good agreement 
as for the amplitude of this d i fferen tial rotation: for the 20 Mo- 
model of iMevnet&Maedeil d2000l) . we find a £Q/Q ~ 0.18 
while their ID-model gives 6Q./H ~ 0.21. 

The dependence of differential rotation with respect to the 
mass of the stars is illustrated in Fig. [14] There we see that the 
trend is that radial differential rotation is decreasing with mass 
while latitudinal variations increase. This is a consequence of the 
variations of the Brunt-Vaisala frequency TV with mass. As mass 
increases, the peak of N 2 near the core decreases and therefore 
the baroclinic torque is weaker, leading to a milder radial differ- 
ential rotation. 

Decreasing metallicity obviously increases differential rota- 
tion as shown by Fig. [15] The stronger latitudinal differential 
rotation is conspicuous. We understand this effect as a conse- 
quence of the reduction of the opacity when metals are removed. 
The gas is then more efficient at heat transport and therefore tem- 
perature gradients are steeper, making the baroclinic torque and 
the differential rotation stronger. 

Finally, we tried to evaluate the effects of time evolution 
along the main sequence. For early-type stars, the variations of 
the hydrogen mass fraction in the convective core can serve as 
a proxy of this evolution. Fig. [16] (top and middle) clearly show 
that at a given fraction of the break-up angular velocity, hydro- 
gen burning leads to an increasing radial differential rotation and 
a decreasing latitudinal variation. This is related to the build- 
ing up of a density jump at the core-envelope interface as dis- 
cussed above. As a consequence, the angular velocity shows a 
rapid variation on this interface as shown in Fig. [161 (bottom). 



5. Conclusions 

In this paper we have presented the first realistic two- 
dimensional models of a rotating star where the differential 
rotation and the meridional circulation are calculated self- 
consistently with the stellar structure. In particular, we have de- 
vised a method in order to compute the velocity field without 
having to compute the Ekman boundary layers. These layers are 
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1.05 
1.00 



Fig. 11. Latitudinal differential rotation at the surface of a 3 M Q 
ZAMS star for several values of the equatorial rotation velocity 
(Q eq ). The group of lower curves shows the latitudinal variations 
of angular velocity at the core-envelope boundary. 



In 



Fig. 12. Latitudinal variations of temperature on an isobar at var- 
ious depth. The change from a decreasing function to an increas- 
ing one occurs near the isobar which crosses the rotation axis at 
a fractional radius of r ~ 0.6. 



indeed extremely thin and would require too much spatial reso- 
lution to be included in the solution. The account of the boundary 
layer is however necessary to lift the indetermination of the in- 
viscid solution with respect to geostrophic flows. The present 
solution therefore includes viscosity a minima so as to make 
the differential rotation well-defined and to insure the balance 
of angular momentum flux between meridional advection and 
diffusion. This simplification does not allow the inclusion of a 
Stewartson layer which requires the full viscosity terms around 
the core. 

The accuracy of the models have been certified through var- 
ious tests: 



- the spectral convergence 

- the monitoring of round-off errors 

- the virial test 





— n«,=o.2fi, 






-■- O. cq =0.5fi, 






n,=o.7n 4 






-- ^=0.9^ 










Fig. 13. Radial profile of differential rotation at the equator (top) 
and the pole (bottom) of a 3 M Q star for several values of the 
equatorial rotation velocity (O e q)- 



0.90 
0.88 










-- n„=o.8n i 



M (M e ) 




M(M 



Fig. 14. Top: Radial differential rotation, taken at equator, for 
ZAMS (homogeneous) models as a function of mass. Bottom 
: surface differential rotation for the same stars. The chemical 
composition is solar. 



- the energy balance test 

The first three tests usually show relative errors less than 
1(T 10 while the energy balance yields much higher errors of or- 
der 10~ 5 because of the rapid variations of opacity. Of course, the 
computation of models close to break-up is less precise because 
of the development of the equatorial cusp (making ^-derivatives 
singular). 

The resulting models have been compared to one- 
dimensional models computed with well-tested codes like 
CESAM or TGEC in the non-rotating case. Slight differences 
of the order or less than a percent are perceptible. They most 
probably come from our use of an analytic formula for the nu- 
clear heating instead of a nuclear network as is standard in ID 
codes. 

Finally, we have compared our models to observational data 
of some rapidly rotating stars: namely a Lyr, a Oph and a 
Leo. These stars have recently been studied in detail with in- 
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Fig. 15. Same as Fig.[l4]but stars with Z=0 metallicity. 




Fig. 16. Top: Variation of the radial differential rotation with 
mass fraction of hydrogen in the core. The value in the enve- 
lope is held fixed at 0.7, while that of the core is decreased as 
would be the case when the stars evolves. The angular velocity 
is half of the break-up: Q. eq = 0.5Qa . Middle: Same as top but for 
the surface differential rotation. Bottom: The ratio of the angular 
velocity at the pole of the core and at the bottom of the envelope 
on the rotation axis. 



terferometers operating in the visible or infrared. Our models 
nicely match the observationally derived parameters of these 
stars suggesting that they are not far from reality and confirm- 
in g the previous good matching of t he gravity darkening found 
bv lEspinosa Lara & Rieutordl d201 lb . 

From these reassuring results we have studied the properties 
of the differential rotation of baroclinic origin that pervade the 
radiative envelope of intermediate-mass stars. The results show 



that it remains small for a homogeneous star: a few percent in 
latitude at the surface and less than 25 % radially. Differential 
rotation seems to decrease with increasing mass and metallic- 
ity. Besides, core evolution increases the radial differential ro- 
tation because of the torque coming from the density jump at 
the core-envelope interface. As shown in Fig. [16] by the end of 
the main sequence, the core rotation may reach three times that 
of the envelope. This result is important in view of the interpre- 
tation of the core rotation of red giants stars descending from 
intermediate-ma ss stars. This rotatio n is now available from as- 
teroseismology dMosser et al. 1 12011 . Before the advent of a full 
account of dynamical evolution in two-dimensional models, our 
results suggest that baroclinic torques may lead to red giants with 
fast rotating cores. 

The foregoing two-dimensional models give a new view of 
the structure and the dynamics of rapidly rotating stars, but 
they raise many new questions as well. We have mentioned 
that of the temperature profile on isobars which behaves non- 
monotonically and for which no simple explanation was found, 
but the most challenging questions are certainly those around 
the core-envelope interface including the Stewartson layer that 
inevitably develops along the tangent cylinder of the convective 
core. Solutions require the account of viscosity and the diffusion 
of elements. However, Ekman layers still need to be avoided 
because they are much too thin to be solved numerically. The 
global solution likely necessitates a combination of the asymp- 
totic solution that we derived in this work and a numerical so- 
lution of the full viscous operator in the interior of the star. 
The density variations of the core-envelope boundary induced 
by core evolution should also be smoothed by some diffusion 
effects; they nevertheless lead to a strong gradient of angular ve- 
locity that emphasizes the mixing of the core and the base of 
the envelope. One may naturally wonder if this dynamical fea- 
ture would lead to some faked core overshooting that is confused 
with a true overshoot. 

Finally, the next challenging issue in this modelling of 
rapidly rotating stars will be the inclusion of time evolution so as 
to monitor the evolution of the distribution of angular momen- 
tum and determine the consequences of the famous rotational 
mixing which is at the heart of the speci fic chemical evolution 
of rotating stars dMaeder & Mevnetll2000t) . 
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Appendix A: Boundary layer equations 

A. 1. Differential operators using boundary layer coordinates 

The boundary layer coordinates (£, t, ip) defined in Sect. I3.4I 
form a set of orthogonal coordinates with scale factors 



h(. = s 



R(t) 



cos 6(t) 



(A.l) 



hu, = R(t) sin t 



It is useful to define a new coordinate s that depends on t 

(A.2) 



only 

s = R(t) sin(r) 
The relation between s- and ^-derivatives is given by 
d cos 5 d 



ds R cos(r - 6) dr 
and the associated scale factor is 
1 



h - 



cos(t - <5) 



(A.3) 



(A.4) 



Based on the scale factors, the common differential operators 
can be expressed as: 
Gradient: 



= -— £ + cos(t-(5)— t+ ~^-<p 
e ot, os s dip 

Divergence: 

_ 1 dv £ cos(r - 5) d , , 1 <9vv 

V • v = - — - + — (sv T ) + - — - ^ - 

s o^ 5 os s dip 



(A.5) 



Curl: 

(cos(t-5) d 1 3v T \- 

s os s dip ) 

' 1 <9vf 1 dvy \ 

K s dip s dij ) 

-77 -COS(T-5)— \<p 
[E Ot, OS J 



Laplacian: 

1 d 2 (p cos(t -5)8 



(A.l) 



e- d£ z s 
Material derivative: 



+ Scos(T _ <5)_L + --^(A.8) 
os \ os ) s 1 dip 1 



i d 2 



(«-V)v= — — +cos(t-S)u t — + — — 

E Ot, OS s Oip 

— ^7+ cos(r - 5) \u T — - + -=- 7- 

E Ot, \ OS s j s dip 

— -f +cos(t-6) \u t — r + +75-7^ 
e ot, \ os s J s dip 



A.2. Vorticity equation 

Let us recall the general form of the vorticity equation 



f (A.9) 



Vx(wx») - s——!p = 

oz 



dQ. 2 A Vp x Vp 



P- 



+ EvAtu + EMJu) (A. 10) 



On the other hand, the interior inviscid solution satisfies the 
equation 



<dQ. 2 \ „ _ VpxVp 
d z /int. 9 P 2 



(A. 11) 



Subtracting the two equations, we have 

Vx(wxu)- 2Qi(z • i)/3ip = EvAco + EM^u) (A. 12) 

where we have used 

Vn = (VQ) mt .+/?f (A. 13) 

Within the boundary layer, the two components (u^ , M T ),of 
the velocity field can be written as 



= B 2 U + £; ■ H;, 
U T = SV + T ■ Hi nt 



(A. 14) 
(A. 15) 



where the interior solution H; nt is ~ 0(s 2 ). The associated vor- 
ticity is 



w = Vx«=w^ = | — + O(e) I <p 



(A. 16) 



It can be shown that the first term in (IA. 12b is 

V x(ioxu) ~ O(s) (A. 17) 

while the last term depends on second derivatives of the velocity 
field, so 



(A.6) EM^u) ~ O 



E\\u\ 



O(e) 



(A. 18) 
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However, as we will see, the remaining term EvAcj is ~ 0(1). 
Indeed 

EAco = EA (cjip) = E |Au - — J <p 

/<9 3 V \\ (A-19) 

-\w +0{s r 

Introducing this in (IA. 1 2b and rearranging terms we obtain, up 
to leading order in e 



d 3 V 



2QS/3 cos(r - 6) + v— = (A.20) 



where we have used the fact that, in the boundary layer, s - s 
and z ■ £ = cos(t - 6). 

A3. Transport of angular momentum 

We will start with the general equation in its compact form 

V ■ (ps 2 Qu) = EV ■ (jUi 2 VQ) (A.21) 

Both terms in the equation are ~ 0(s 2 ) for the interior variables, 
but they are ~ 0(e) within the boundary layer. Let us start with 
the right hand side 

EV ■ (jus 2 VQ.) = e 2 V • (// <r(VQ) int ) + e 2 V ■ (jjs 2 /3$) 

= e^s 2 P ) + 0(e 2 ) ^ 

Note that the term V ■ (fis 2 ( VQ)i nt .) depends only on interior vari- 
ables and so it is ~ 0(1). The vertical derivative of a interior 
variable with respect to £ is ~ 0(e), so we can write 

£V-(^s 2 Vn)=s^ 2 f +0(s 2 ) (A.23) 

Now, we calculate the remaining term 

V ■ (ps 2 Qu) = ps 2 u ■ VQ + s 2 QV • (pu) + 2sQps ■ u (A.24) 
From the continuity equation, we know that V • (pu) = 0, then 

V ■ (ps 2 Q.u) = ps 2 u ■ VQ + 2sQ.ps ■ u (A.25) 

Substituting VQ = (VQ) int + jfif and u = u mi + s 2 U% + sVf, we 
have (remember that u mt ~ 0(e 2 )) 

V ■ (ps 2 Qu) = eps 2 Vf ■ (VQ) int + 2ssQ.p(s ■ f)V + 0(s 2 ) (A.26) 

Substituting dA.26t and (IA.231 > in dA.2U , up to leading order in 
s, we get 

sVf ■ (VQ) int + 2Q(S ■ f)V = vs% (A.21) 

We know that, in the boundary layer s — s and 

f ■ (VQ) int = cos(t - 6)-£ (A.28) 
os 

s • f = cos(r - (S) (A.29) 
then, rearranging terms we obtain the final form 

hn + s^)vcos(T-S)-sv^ = 0. (A.30) 
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